ON ULTRASONIC SIGNAL PROCESSING AND ITS 
APPLICATION TO TOMOGRAPHIC RECONSTRUCTION 


by 

KANDREGULA PRASAD 


TH 

netem/n 

Ps«e^ 



department of mechanical engineering 
Indian Institute of Technology, Kanpur 

April, 2001 



ON ULTRASONIC SIGNAL PROCESSING AND ITS 
APPLICATION TO TOMOGRAPHIC RECONSTRUCTION 


mm 


A Thesis Submitted 

in Partial Fulfillment of the Requirements 
for the Degree of 

MASTER OF TECHNOLOGY 


by 

KANDREGULA PRASAD 



to the 

DEPARTMENT OF MECHANICAL ENGINEERING 
INDIAN INSTITUTE OF TECHNOLOGY KANPUR 

APRIL, 2001 



2 o JUL 200l / h B 

g w t f re ^Tfr?TT5T 
Wcffr sfaflffaft tfWT 
dferifer ^o A..134t1i> a^®«La® ®S^ 



CERTIFICATE 


It is certified that the work contained in the thesis entitled ON ULTRASONIC 
SIGNAL PROCESSING AND ITS APPLICATION TO TOMOGRAPHIC 
RECONSTRUCTION by Kandregula Prasad has been carried out under my 
supervision and that this work has not been submitted elsewhere for a degree. 


April, 2001 


|\i aj 

Dr. N. N. Kishore 
(Professor) 

Department of Mechanical Engineering, 
Indian Institute of Technology, Kanpur. 


ii 



DEDICATED 

TO 

MY BELOVED PARENTS 


ACKNOWLEDGEMENTS 


I wish to express my profound gratitude and indebtedness towards my thesis supervisor 
Dr. N.N. Kishore for his inspiring guidance, invaluable suggestions and constant 
encouragement. 

I wish to extremely thankful to Mr. S. K. Rathore for his valuable suggestions 
and constant encouragement towards completion of my Thesis. 

I wish to express my thanks to Dr. V. Raghuram for his suggestions and 
encouragements. I appreciate and extend my thanks to my lab mates Dr. A. K. 
Agrawal, Shashi bhushan, Samuel, Bnrao, Dali raju, Sarin and Trivedi. 

I would like to thank all of my friends Swaroop, Kovvila, PLN, Ganesh, 
Narendra, Ramesh, Lachi, KNS, Challa, Bhanu, Anil and Laxman for making my stay 
at IITK very enjoyable and memorable. I will cherish the moments forever. 

I wish to record my special thanks to my uncle Mr. P. Ramana babu and my 
cousin for motivating me to pursue my M. Tech program. I wish to extend my thanks to 
my sisters and cousin. I am also very much thankful to my Uncle Mr. M. Apparao for 
his valuable help. 


Indian Institute of Technology, Kanpur. 
April 2001 


Kandregula Prasad 



ABSTRACT 


In Ultrasonic Non-destructive testing and Evaluation (NDE), the different features of 
ultrasonic waves, namely frequency, amplitude, time of flight etc. can be used for 
detection and sizing of internal defects. Ultrasonic Tomography (UT) is used to 
characterize a defect or in-homogeneity by Time of Flight (TOF) information. The 
accuracy of reconstruction depends on the correctness of TOF determination. The present 
work is on the use of signal processing techniques for TOF measurement and consequent 
image reconstruction. In practice, the ultrasonic beams in an elastic medium is generated 
by one probe (transmitter) and received through another probe (receiver). An Automated 
Ultrasonic Experimental Setup was used for data acquisition in an A-scan format. 
Discrete Wavelet Transform (DWT) was used for noise suppression from all the acquired 
signals, by using ‘db5’ mother wavelet. Samples were prepared using combinations of 
Perspex and Nylon as base materials and Epoxy resin and Aluminium as insert materials. 
When a pulse wave is received through an interface, the signal has overlapped modes. 
The modes of these overlapped (i.e. multi-path) signals were separated by using least- 
squares method and cross correlation method. Criterion for selection of a mode was based 
on the number of reference signals available. Based on these results, an energy-based 
criterion has been proposed for the selection of desired mode from the separated modes. 
The Time of Flight of these modes have been used for tomographic reconstruction of a 
specimen having defect using a Multiplicative Algebraic Reconstruction Technique 
(MART) algorithm. The reconstructed images for 2- view and 4-view projections are 
presented and it has been shown that for reconstruction of true shape, 2-view-projections 
are not sufficient. The effect of relaxation parameter (X) on the reconstruction accuracy 
and convergence rate has also been studied. 
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Chapter 1 
INTRODUCTION 


1.1 INTRODUCTION 

Non-Destructive testing and Evaluation (NDE) techniques play a major role in quality 
assurance of structural members during manufacturing stage and operating life. Being 
one of the most commonly used Non-destructive Testing (NDT) methods. Ultrasonic 
Testing (UT) is fast developing in the recent years. The UT method uses ultrasonic waves 
for detection and sizing of internal defects of materials. In this, piezo-electric transducer 
(probe) generates ultrasonic waves, which propagate in the elastic medium and are 
detected either by the same or by a different transducer. 

In general the transit time of the acoustic echo wave is used to measure the 
distance of a flaw/insertion from the probe and the amplitude of the reflected wave is 
used to size the flaw. In the pulse echo technique, a single probe is used to transmit and 
receive the signal. Another technique called through transmission employs two separate 
transducers (a transmitter and a receiver) to transmit and receive the wave. This technique 
requires access to both sides of the specimen, as the probes are placed on opposite sides 
of the specimen. 

Ultrasonic signals can be displayed in A-scan, B-scan or C-scan formats. An A-scan 
indicates a variation of signal amplitude with time at a point. The familiar screen image is 
on A-scan presentation, which furnishes one-dimensional description of a given test 
point. In the case of B-scan, the displayed signal is a time versus linear position trace. In 
case of two-dimensional (area) scanning of test piece, the test results can be represented 
by means of a C-scan. This furnishes a top view of the scanned surface. C-scan systems 
are usually computer controlled and the image displayed is colour coded to facilitate 
interpretation. Ultrasonic systems displaying both B-scan and C-scan formats are called 
P-scan systems. 

Wavelet analysis in ultrasonic nondestructive applications has been shown to be a 
useful tool for the interpretation and enhancement of ultrasonic data in the context of 
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nondestructive evaluation. The Wavelet Transform (WT) divides the signal into different 
frequencies/scales/levels in terms of basis functions obtained by compression/dilation and 
translation of mother wavelet. Compared to other time-frequency representations like 
short time Fourier transform and Wigner-ville transform, the wavelet transform provides 
spectral representation and temporal order of the signal components simultaneously. 
Another important property is that the signal reconstruction does not involve global 
averaging in time or frequency domains because of good localization of the wavelet 
coefficients in both domains. 

In ultrasonic techniques, information on defect characterization requires more 
evolved techniques than classical methods. Modern non-destructive testing of materials 
has to provide the highest possible detection probability, the correct size and exact 
orientation of defects in the specimen. The simplest mechanical acoustical 
characterizations of materials deal with velocity measurements. Due to mode conversion 
of the solid-solid/ liquid-solid interface, two quasi waves coexist in the solid medium. 
The spatial and temporal separations of the two modes are limited in certain conditions. 

The technique of Computerized Tomography (CT) has established itself as a leading 
tool in diagnostic radiology over past thirty years, and is being applied in the NDE of 
engineering structures in a variety of situations. Image reconstruction from projections is 
based either on transform methods or on series expansion methods. Series expansion 
methods are iterative. 

This study is on the reconstruction of the samples with artificially implanted defects 
is carried out using Ultrasonic Computed Tomography (UCT). When a pulse wave is 
received through a region of in-homogeneity interface, the signal may have two or more 
overlapped modes. The accurate and correct TOF data estimation from these overlapped 
signals demands the application of appropriate signal analysis techniques, which may 
lead to good reconstruction quality. 

1.2 LITERATURE SURVEY 

In recent years, many mathematical tools have been developed for the analysis of data in 
Non-destructive evaluation. Wavelet analysis is one such important technique. Roy et al 
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[1] have studied the continuous wavelet transform (CWT) applications for signal analysis 
wherein a disturbing echo of high amplitude is extracted from an experimental ultrasonic 
signal, which is resulting from a clad state of surface or a side-drilled hole of steel block. 
Agostino et al [2] utilized wavelet transform for noise suppression to improve ultrasonic 
flaw detection. The improvement in detection was experimentally verified using steel 
samples with simulated flaws. Guang-Ming et al [3] performed computer simulation to 
verify the signal detection improvements for an ultrasonic wave embodied in white noise 
using cast steel samples with artificial flaws. Later a theoretical model to estimate the 
optimal frequency-to-bandwidth of the Gaussian wavelet for denoising in ultrasonic 
NDE is presented [4] with experimental verification. 

Castagnede et al [5] analyzed the optimum possibilities offered by correlation 
method for time of flight measurement when mixing of different modes is present. A 
criterion for the separation of the different modes, the systematical error associated with 
their overlapping and an algorithm for unraveling them are also described. Drai et al [6] 
have developed some signal processing algorithms to facilitate real time processing of 
ultrasonic NDT data. These algorithms include applications of cross correlation function 
for planar and volumetric defect discrimination. 

Mueller et al [7] discussed several methods of ultrasonic tomography with 
diffraction effects. Kak [8] reviewed the major developments in imaging with computed 
tomography (CT) using X-ray emission and ultrasound sources. Errors in reconstruction 
are used to compare different algorithms, are used by Subbarao et al [9], and it is also 
found that noise level in the projection data increases the error and slows down the 
convergence rate in the reconstructed field. The maximum error and RMS error have 
been found highest in ART and lowest in the MART for given projection data. 

1.3 PRESENT WORK 

In the present work, the experimental data obtained from an ultrasonic automated 
scanning setup has been used for tomographic reconstruction of defects of different 
shapes in various samples. One dimensional discrete wavelet transform was applied for 
noise suppression from the acquired waveforms. The waveforms received through 
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inhomogeneity contain two overlapped modes. For separation of these overlapped modes, 
least-squares method and cross correlation method are presented. A measure of energy- 
based criterion is suggested to identify the desired mode from the separated modes, 
whose TOF value is to be considered. This processed time of flight (TOF) is used to 
reconstruction of the image using computed tomographic (CT) algorithms. 

In chapter 2, the details of the experimental setup and data acquisition procedure 
are discussed. 

Chapter 3 presents the wavelet analysis and signal analysis techniques. In this 
Chapter wavelet application for noise suppression is also given. 

Chapter 4 deals with mode separation of multi-path signal using least-squares 
method and cross correlation method. Time of flight estimation with various criteria is 
also given in this chapter. 

Chapter 5 presents theories of various algebraic reconstruction techniques of 
computer tomography. Reconstruction of the various domains of the tested samples is 
also discussed with their images in terms of slowness. 

Chapter 6 gives the conclusions of the present work and suggestions for the future - 


work. 



Chanter 2 

EXPERIMENTAL SET-UP AND PROCEDURE 

This Chapter briefly describes the details of the experimental set-up and the aspects of 
data collection. 

2.1 EXPERIMENTAL SETUP 

The schematic diagram of ultrasonic automated scanning setup is shown Fig.2.1, to 
perform A-scan. The present ultrasonic experimental setup, developed in-house 
(Mahajan [26], Venkatesh [27]), consists of scanning tank, a high resolution type 
Ultrasonic Flaw Detector (UFD), a high speed data acquisition and other accessories 
for precision movements of the transducers, data collection and storage. A brief 
description about the setup is given in the following paragraphs. 

2.1.1 SCANNING TANK 

The purpose of the scanning tank is to house the stepper motors, the probes and 
specimen to be tested. As is common practice, the probes and specimens are 
immersed in water, which acts as a coupling medium. The tank is made up of Perspex 
plate. Two lead screws, each having pitch of 4 mm, have been fixed in mutually 
perpendicular directions on a frame located at the top of the tank. The probes mounted 
on the lead screw could be moved in x-y plane with a step of 1 mm in either direction. 

2.1.2 ULTRASONIC FLAW DETECTOR 

The Ultrasonic Flaw Detector (UFD) (Krautkrammer-Branson make, USIP12) is used 
for generating and collecting the ultrasonic signal with the help of probes. The UFD 
can be used with a single probe in pulse-echo transmission method or with two probes 
in through-transmission method. The received signal can be obtained from UFD in 
rectified or radio frequency (RF) mode. There are four controls (for Gain adjustment) 
by means of which the gain of the received signal can be adjusted to the desired level 
to reduce noise and avoid any saturation. Various controls are provided for fine 
adjustments of frequency range, echo shape, pulse repetition frequency, Distance 
Amplitude Correction (DAC), Monitor adjustment, back wall echo. 



2.1.3 COMPUTER 

A microcomputer (PC- AT 386SX) is used to control the stepper motors, UFD, A/D 
converter for data collection and storage. 

2.1.4 STEPPER MOTOR AND CONTROLLER 

The stepper motor (STM 601 with 12v input) of torque capacity 2kg-cm with 
maximum step rate of 1023 steps/sec were used. The motor can be moved either 
clock- wise or anti-clock wise direction through stepper motor control card PCL-211 
(manufactured by Dynalog Microsystems). In present setup motors are moved with 50 
clock pulses at a time, which corresponds to 1 mm linear movement of the probe. 

2.1.5 ANALOG TO DIGITAL CONVERTER 

SONTEK STR* 8100 A/D board is used to digitize RF signal from UFD with a 
maximum data acquisition flexibility through a 64k on-board high-speed memory 
buffer. The various board functions, under software control, are input channel 
selection AC/DC coupling, input voltage range, RF/Video made, sampling rates, 
trigger selection, clock control, threshold phase and level, board selection and 
interrupt enabling. The board’s high-speed data acquisition memory is mapped 
directly into PC memory space. Once the waveform has been captured, the PC can 
transfer data off the board at 1.5 MHz rate (8 bit) or 3 MHz rate (16 bit) using simple 
memory block instructions. This board provides freedom from the static architecture 
of a stand alone instrument and the computer. The D100A Digital Oscilloscope is 
used to examine the waveform for selecting the appropriate scan parameters such as 
delay, gate start, gate length etc. In the present setup the A/D board is installed outside 
the PC using the PCX-795 PC XT/AT bus expansion system. This board has on-board 
1,2,4,8,16 and 32 wait state generation, switch setting for memory bank selection, 
switch setting for I/O selection and switch setting for DMA channel selection. 

2.1.6 TRANSDUCERS 

The transducers (probes) used in the present investigation are Immersion Straight 
Beam probes Z10M, which generate broadband longitudinal waves. They are 


connected to the UFD by waterproof connectors. The fixtures were fabricated to hold 
the transducers in vertical position. 

Immersion probes are mainly applied for mechanized or automatic ultrasonic test. 
In principle they work same as contact probes. A major part of the tests is made in 
immersion tests tanks filled with water. The diameter of the transducer element is 5 
mm (element size). The size of the element strongly affects the shape of the 
transmitted sound field. The probe has nominal frequency of 10MHz. The frequency 
has a great influence on the evaluation of reflectors with increasing frequency, the 
echo height from non-vertically positioned reflectors to the sound beam decreases. 


2.2 EXPERIMENTAL PROCEDURE 

In the present work, the various specimens tested are given in the following Table. 



Material 


Table 2.1 


Longitudinal wave Velocity (m/s) 


Body 


Perspex 


Ny 


Perspex 


Nylon 


Perspex 


Insert 


Epoxy 


Body 


Insert 


Fig. No. 


2.2(a) 




Perspex 


Epoxy 


Aluminium 


Epoxy 


Epoxy 


2.2(c) 


2.2(d) 


2.2(e) 


2.2(f) 


* Epoxy resin is LY556 epoxy cured with 10% by weight of hardener HY951, 
supplied by Ciba-Geigy of India ltd. 

The specimen to be investigated is cleaned and area under investigation is 
marked. Then, the specimen is positioned between the Transmitter and Receiver as 
shown in fig.2.3. Probes are brought to start position by moving the motor. The gap 
between the specimen plate surface and the probe surface are adjusted, so as to 
distinctly observe the initial burst on the UFD oscilloscope screen. The gain 
adjustment controls are so adjusted that the peak of the signal does not exceed 80% of 
the maximum screen height of the UFD to avoid saturation. The RF signals from UFD 























can be also viewed on the digital oscilloscope provided by the software of the A/D 
board. The gate parameters on the signal to be digitized are made in data file. In these 
experiments, the data collection is in one dimension only. 

A master control code (in C-language) has been developed for the present 
investigation controls precise movements of stepper motors and controls the A/D card 
to collect received waveform in the digitized form. The code is also capable to extract 
different time and frequency domain features of the waveform during scanning. The 
initial settings, such as channels, AC/DC coupling, clock source (Intem/Extem), 
sampling rate (100/200/300/400), trigger type (off/on), post trigger delay, buffer start, 
buffer length, horizontal delay, gate start, gate length, xfer type (8/16) are set in a data 
file prior to the execution of the code. 

2.2 CLOSURE 

In this Chapter the details of experimental setup and procedure were presented. 
Different components like scanning tank, ultrasonic flaw detector, computer, stepper 
motor, stepper motor controller, analog to digital converter and transducer were 
explained. Experimental procedure for data acquisition was also explained. 



T: Trans mitter Probe 



Fig. 2.1. Schematic Diagram of Ultrasonic Automated Scanning Set-up 
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(a) SAMPLE 1 
A— Epoxy Resin, B— Perspex 
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(b) SAMPLE 2 
A— Epoxy Resin, B— Nylon 
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(c) SAMPLE 3 (d) SAMPLE 4 


A— Epoxy Resin, B— Perspex A— Epoxy Resin, B— Nylon 

Figs. 2.2(a-d). Different Type of Samples used for the Present work 
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(*) SAMPLE 5 
A- Aluminium, B- Perspex 


(f) SAMPLE 6 
A— Epoxy Rcnn, B— Nylon 




(g) SAMPLE7 
A— Epoxy Resin, B— Perspex 



(b) Setup tot collecting DATA at 45 and 135 
degrees foe SAMPLE 7 
A— Epoxy Resin, B— Perspex 


Figs. 2.2(e-h). Different Type of Samples used for the Present work 
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Fig. 2.3. Positions of Transmitter and Receiver for data acquisition 
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Chapter 3 

NOISE SUPPRESSION USING WAVELET TRANSFORM: 

A SIGNAL PROCESSOR 

3.1 INTRODUCTION 

Signal processing has become an essential part of contemporary scientific and 
technological investigation. Signal processing is used in telecommunications, in the 
transmission and analysis of satellite images, and in medical imaging (echo graph, 
topography and magnetic resonance), all of which involve the analysis and interpretation 
of complex time series. In the recent times, this technique is also being applied in Non- 
destructive Evaluation (NDE) of materials. 

The specific goals of signal processing in NDE are. 

• To improve inspection reliability, 

• To improve defect detection, 

• To improve defect characterization, and 

• To generate information about the material properties to asses residual life of a 
structure. 

The Wavelet Transform (WT) is a recent mathematical technique for processing signals 
with time-varying spectra, has excellent characteristics of time-frequency localization. 
Wavelet analysis has wide applications in earthquake prediction, signal processing, image 
compression, turbulence, non-destructive evaluation (NDE) etc.. In the following 
sections, the basic operations of wavelet transform are presented. 

3.2 BASIC SIGNAL ANALYSIS TECHNIQUES 

3.2.1 FOURIER ANALYSIS 

Fourier analysis is well known technique of signal processing, which breaks down the 
signal into constituent sinusoids of different frequencies. 
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The Fourier Transform (FT) is a technique to transform time-domain signal f(t) 
into frequency- domain signal F(co). Fig.3.1 shows an example of the Fourier Transform. 
Mathematically the Fourier Transform is given by 

n<*>)=j/( t).e-.dt (3.1) 

—00 

and its inverse, 

oo 

f(t)= \F(co).e m .dco (3.2) 

—oo 

where, oo is the frequency. 

In Fourier Transform, the main drawback is that the time information is hidden. 
When looking at a FT of a signal, it cannot be seen when the particular spectral 
components appear exactly in the signal. 

3.2.2 SHORT TIME FOURIER ANALYSIS (STFA) 

Short-Time Fourier Transform (STFT) is a modification of the Fourier transform. STFT 
has adapted by Dennis Gabor (1946), to analyze a small section of the signal at a time, a 
technique is called widowing a signal. Fig. 3.2(a) shows transformation of a signal into 
two-dimensional function of time and frequency. Mathematically the STFT, with window 
function w(t-r) is expressed as 

oo 

STFT (oo) = \f{t).w{t~T).e ,M .dt (3.3) 

-OO 

where, w (t-r) is the analysis filter or analysis window. 

The STFT represents a compromise between the time-domain and frequency- 
domain views of a signal. It provides information about frequency content in the window 
of the signal at the given time. However, this information is obtained with limited 
precision, and that precision is determined by the size of the window. The drawback is 
that once a particular size is chosen for the time window, that window is the same for all 
frequencies. Many signals require a more flexible approach, where the window size can 
vary to determine more accurately in either time-domain or frequency-domain. 
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3.2.3 WAVELET ANALYSIS (WA) 

Wavelet Analysis is the next logical step in which the windowing technique with variable 
sized regions is used. Wavelet analysis allows both low frequency information with long 
time intervals and high frequency information with shorter regions. Fig. 3.2(b) shows that 
wavelet analysis does not use a time-frequency region, but rather a time-scale region. 
Mathematically, the Wavelet Transform (WT) of a signal is defined as follows: 


W f (a,b) = J/(<) 



t -b 
a 


L dt 


(3.4) 


where, \g(t) is a mother wavelet, b is a time shift parameter and a is a scaling parameter. 
By definition, the WT is the correlation between the signal and set of basic Wavelets. 

Wavelet analysis produces a time-scale view of a signal, with scaling and shifting. 
Scaling a Wavelet means stretching (or compressing) it. Fig. 3.3 shows wavelets with 
different scaling parameters. Shifting a Wavelet simply means delaying (or hastening) its 
onset. Mathematically, delaying a function f(t) by t is represented by f(t-z). Fig. 3.4 
shows wavelets with different shift factors. 

Wavelet analysis is capable of revealing aspects of data that other signal analysis 
techniques miss such as trends, break down points, discontinuities in the higher 
derivatives and self-similarity. Wavelet analysis can compress or de-noise a signal 
without appreciable degradation. 

Wavelet analysis can be performed in two ways 

1 . Using continuously translated and dilated versions of mother wavelet (CWT), or 

2. Using discretely translated and dilated versions of mother wavelet (DWT). 


3.2.3. 1 Continuous Wavelet Transform (CWT) 

The continuous wavelet transform (CWT) is defined as the sum over all time of the signal 
multiplied by scaled and shifted versions of the wavelet function % which can be 
represented as 
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(3.5) 


C(a,b)= \f(t)y/' a b (t) dt 

-CO 

where, C (a, b) are the coefficients of the wavelet transform, f (t) is the signal being 
analyzed and * denotes the complex conjugate of y / ab . Here, y/ a b is the shifted and 

scaled version of the mother wavelet y/(t) as given by, 

— |, ( a,b)eR 2 , a > 0 (3.5a) 

bj 

The word continuous in CWT means, that the transform operates the wavelets on 
the signal with continuous in scaling and shifting. During the computation, the analyzing 
wavelet is shifted smoothly over the full domain of the analyzed signal/function. 



3.23.2 Discrete Wavelet Transform (DWT) 

Calculating W avelet coefficients in CWT at every possible scale generates a very large 
amount of data. This disadvantage of CWT is overcome by choosing scales and positions 
based on powers of two, so called dyadic scales and position, which is more efficient and 
just as accurate. Such an analysis is called Discrete Wavelet Transform (DWT). The 
generation of the wavelets and calculation of the DWT are well matched to digital 
computer. 

The wavelet transform can also be expressed in discrete form 


where, 


(3.6) 

m n \ ^0 J 

C m ,n = < n \m ¥ { t ~ na m ° T \dt (3.7) 

V a 0 J 

a = < (3.8) 

b = na™T (3.9) 


where, T is the sampling period, A and ao are constants, where ao=2 for the dyadic 
wavelet transform, m is depend upon the number of decomposition levels and n is depend 
upon the length of the signal to be decomposed. 
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Multi Resolution Analysis (MRA) 

Multi resolution wavelet analysis allows the decomposition of a function/signal in 
progression of successive approximation and details, corresponding to different scales. 
The approximations are high-scale, low frequency components and the details are the 
low-scale, high frequency components of the signal. The difference between the actual 
signal and its approximation of order n is called its residual. Intuitively, the 
approximation is relatively smooth, and detail being composed of high frequency 
components. The detail corresponds to the difference between two successive levels of 
approximation. 

Multiple-Level Decomposition 

In DWT, the decomposition process can be iterated, with successive approximations 
being decomposed in turn, so that one signal can be broken into many lower-resolution 
components. This is called the Wavelet decomposition tree. Fig. 3.5 shows the Multi- 
level (three level) wavelet decomposition tree. 

Wavelet Reconstruction 

Wavelet reconstruction, or synthesis is the process of reconstructing the decomposed 
signal using a set of decomposed coefficients. The mathematical manipulation that effects 
synthesis is called the inverse discrete wctvelet transform (IDWT). 

3.2.3.3 Properties of Wavelets 

The best-performing mother wavelet should have the following three properties. 

1 . The admissible condition is given by 

oo 

jV(/) dt = 0 (3.10) 

—oo 

2. The wavelet has the characteristic of good time-frequency localization. This 
property makes the wavelet transform is powerful tool for analysis of non- 
stationary signals such as sound, seismic, electromagnetic signals etc. 
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3. The power spectrum of the wavelet is well matched with that of the signal to be 
analyzed. 

3.23.4 Different types of Wavelets 

A wavelet is a waveform of limited duration, which has the above properties. Different 
types of mother wavelets are discussed in the following paragraphs. 

The simplest wavelet is Haar wavelet, which is a step function (Figs. 3.6(a) and 
3.6(b)). 

Ingrid Daubechies, one of the brightest stars in the world of wavelet research, 
invented what are called compactly supported orthonormal wavelets. The names of the 
Daubechies family wavelets are written as dbN, where N is the order, and db is the 
“surname” of the wavelet. The dbl wavelet is same as Haar. Figs. 3.6(c) through 3.6(p) 
shows the scaling and wavelet functions of the various Daubechies wavelets. 

In Biorthogonal wavelets two types of wavelet functions are used: one for 
decomposition and other for reconstruction instead of the same one. Coiflets, Symlets, 
Morlet, Mexican Hat, Mayer etc. are other different type of wavelets suitable for various 
purposes. 

There are different wavelets that can be used to decompose the signal and extract 
feature vector. These wavelets include Daubechies (dbl-dblO), Biorthogonal (biorl.3- 
bior6.8), symlets (sym2-sym8), coiflets (coifl-coif5), Morlet, Mexican hat and Meyer. 
The criterion for selection of the proper mother wavelet is based on the shape of the 
wavelet function, which could match well with the shape of the signal mode to be 
analyzed. In the present work, Daubechies, order 5 (‘db5’) has been chosen, because it 
seems to match with the signal mode to be analyzed and also it posses a good time- 
frequency localization. Figs. 3.6(a) to 3.6(p) shows different type of wavelet and scaling 
functions. Fig. 3.6(j) shows the Daubechies-5 (‘db5’) that has been used in the present 
work. 

3.3 NOISE SUPPRESSION 

In Ultrasonic Nondestructive Evaluation (NDE), the observed ultrasonic signal S(t) can 
be expressed as the sum of the two components. 

^W = t(0+«(0 (3-11) 
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where, y(t) is the reflected signal of the defect and n(t) is the noise. Removal noise is 
extremely important in the ultrasonic defect detection to identify the defects correctly. 

The main steps of the signal noise removal are. 

1 . Decomposition of the signal into coefficients, 

2. Separation of the approximation and detail coefficients from the coefficients, and 

3. Reconstruction of the de-noised signal from the approximation coefficients and their 
different levels of the noises from the detail coefficients. 

A code is developed for the noise suppression from the signal is given in the appendix A. 
In this, the analysis is limited to 3-level decomposition of the signal. In 3-level 
decomposition, the function wavedec (in the MATLAB toolbox) returns the coefficients 
of all the components (i.e. third-level approximation and first three levels of details) and 
concatenated into one vector, C. Vector L gives the lengths of each component. In the 
second step, the coefficients of approximation and details have been extracted from 
decomposed vector, C. In the third step, the de-noised signal (approximation) and the 
three levels of noises (details) have been reconstructed from their corresponding 
coefficients. Fig. 3.7 shows various levels of suppressed noises (details) and de-noised 
signal (approximation) of a non-overlapped signal, which is received through an Epoxy 
inserted Perspex sample with 3-level decomposition. As the level of decomposition 
increases, the frequency of the detail (noise) decreases. Figs. 3.8 and 3.9 show the noise 
suppression of the overlapped signals, which are received through region of in- 
homogeneity of samples Epoxy inserted Nylon sample and Epoxy inserted Perspex 
sample respectively. In fig. 3.10, the original signal is an overlap of two modes with very 
small time shift, which is received through a region of in-homogeneity of Perspex sample 
with Epoxy resin rectangular insertion of 20mmxl0mm. In this figure the de-noised 
signal (approximation) looks like a single mode signal, just after 1 -level decomposition. 
■The signal has lost some important information (i.e. two peaks of two different modes). 
This is the limitation for the noise suppression. So, the de-noising characteristics depend 
on the level of decomposition and signal to be analyzed. Therefore, in the present work, 
different levels have been chosen for different signals according to their shapes. 
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3.4 CLOSURE 

In this Chapter, the basis of FT, STFT and Wavelet Transform (WT) have been 
presented. Some of the important mother wavelets and their properties have been 
discussed. Noise removal using multi-level DWT, which has been used in the present 
work, has also been explained. 
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Fig. 3.1 Fourier Transform for transforming a time-domain signal to frequency-domain 




(a) STFT (Gabor) (b) Wavelet Analysis 


Figs. 3.2(a-b) Short-Time Transform (Time-frequency domain) and Wavelet Transform 
(Time-scale domain) 
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Fig. 3.3 Scaling (stretching or compression) of the wavelets with different scale factors 
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Fig. 3.4 Shifting (delaying or hastening) of the wavelets with different shift parameters 
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Fig 3.5 Multiple-Level Wavelet decomposition tree with three-level decomposition 
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(a) dbl/Haar Scaling function (b) dbl/Haar Wavelet function 



(c) db2 Scaling function (d) db2 Wavelet function 



(e) db3 Scaling function (f) db3 Wavelet function 



(g) db4 Scaling function (h) db4 Wavelet function 


Figs. 3,6(a-h) Different wavelets (Daubechies) scaling and wavelet functions 
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Figs. 3.6(i-p) Different wavelets (Daubechies) scaling and wavelet functions 
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Fig. 3.7 


Fig. 3.8 
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Fig. 3.9 Wavelet 3-level decomposition of overlapped signal at border of rectangular 
insertion of 20mmx20mm (Epoxy resin) in Perspex sample with c db5’ 
wavelet 
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Fig. 3.10 Wavelet 1 -level decomposition of overlapped signal at border of rectangular 
insertion of 20mmxl0mm (Epoxy resin) in Perspex sample with ‘db5’ 
wavelet 
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Qha^tejj 4 

RESOLUTION OF INDIVIDUAL MODES IN A MULTI- 
PATH SIGNAL AND TIME OF FLIGHT ESIMATION 

4.1 INTRODUCTION 

In many applications such as radar and sonar data processing, geological acoustic 
sounding, ultrasound based non-destructive testing, and medical imaging procedures the 
determination of certain event times is faced with the problem of resolving unknown 
number of close spaced, overlapping and noisy echoes of a signal with known shape. 
When a signal is passing through solid-solid interface, the signal path splits into two or 
more paths and arrives at the receiver at different times. Thus, the signal is collected at 
the receiver is superposition of these individual signal modes and is called overlapped or 
multi-path signal. In figs. 4.1(a), 4.2(a), 4.3(a) and 4.4(a) show the multi-paths of various 
samples. In all these cases, when the signal is passes through the border of in- 
homogeneity the signal path splits in to two separate paths and reach the receiver and plot 
as an overlapped or multi-path signal with some time-shift. These types of multi-path 
signals have been observed in the present work. Then the received multi-path S(t) signal 
can be described mathematically as 

S{t) = Yj A ' R ( t ~ T >) +n ( t ) (4.1) 

/=I 

where, A, and x,- are the amplitudes and arrival times respectively, for the z'th multi-path 
signal, N is the number of multi-paths, R(t) is representative reference signal from a path 
and n(t) is the additive noise. 

In the previous chapter, about suppression of additive noise has been discussed. 
After the removal of the noise, the de-noised signal S(t) becomes, 

S(0 = XA,R(t-T l ) (4.2) 

1=1 

The experimental de-noised signals at various potions are shown in figs. 4.1 to 
4.4. Fig. 4.1(a) shows a Perspex sample of 50mmx50mm size with square insertion 
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(Epoxy resin) of 20mmx20mm. Fig. 4.1(b) shows received signals at three different 
positions of location shown in fig. 4.1(a). It can be seen that the signal has only one mode 
when it passes through the complete homogeneous material. If it passes through the 
interface region of the in-homogeneity, the signal is a multi-path signal with two signal 
modes as shown in fig. 4.1(b). If the sound velocity through insert material is faster than 
the base material, then the signals are as shown in figs. 4.1(c) and 4.1(d). In fig 4.1(d), the 
two signal modes with more time shift, because the sound velocity difference between the 
base material (Perspex) and inserted material (Aluminium) is very high. When the 
insertion thickness is small the multi path signal are overlapped closely as shown fig. 
4.2(b) for sample shown in fig. 4.2(a). 

In real practice perfect rectangular/square insertions exist rarely. So, the present 
work focuses on elliptical insertion as shown in fig. 4.3(a). Figs. 4.3(b) and 4.3(c) show 
signals for various positions of the transmitter for the case of Epoxy resin inserted in 
Nylon sample and Perspex sample with Epoxy resin insertion. 

Time of flight (TOF) is the important parameter in NDT. In the present work, for 
separating the different modes of the multi-path signal, two different types of algorithms 
have been used. The methods are 

1 . A Least-squares method algorithm for resolution of multi path signal, and 

2. Cross correlation method using Fast Fourier Transform (FFT) has been used, when 
only one reference (single path) signal is available. 

4.2 LEAST-SQUARES (LS) METHOD 

The method of Least squares is systematic procedure to fit a unique curve through given 
data points. In the present work, the multi-path signal (overlap of two modes) is separated 
with the help of reference (single path) signals. 

Let the received signal S(t) be assumed to consist of two modes. Then S(t) can be 
represented as, 

S(t) = A,R l (t)+A 2 R 2 {t)+A 3 (4.3) 
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where, Rj(t) and R 2 O) are two reference signals, A; and A 2 are amplitudes of the modes, 
and A 3 is DC-offset of the signal from the time axis. If e } is the error of approximation at 
j th point, then 

e,=S,-(4(i, l+AM, + 4 ) (4.4) 

If the signal is descritized at M points, then total square error equation can be written as, 

£ = Ek -UwW), +4)] 

7=1 


= E4 (4-5) 

7=1 

Then, the method of least squares consists of determining such that the error E is 
minimized. 

For E to be minimum, the conditions are 

fr-o-f :- 2 ( R > l I s j - U fa \ + 4 fa ), + 4 1 . («•«■») 

OA\ ; .i 

= 0 = f>2(7! 2 ),[s, -(4(7!,), +4(7!,), + 4 )], and (4.6b) 

7=1 


^- = 0 = E-2fe -(4(7!,), +4 (7!,), +4l 

The above equations in the matrix notation, is given as 

MM-W 


(4.6c) 


(4.7) 


where, 


M= 

M1M2M3 

F F F 

r 2\ r 22 i 23 

. w- 

to *-* 

► , and {g} = < 

'N 

o o 


F F F 

L/ 31 M2 M3 J 


A 

V - 1 J 


1^3 J 


where, 

M IA M 

fu = Efa )) , Fn = Efa)f F » = El = JW • 

<=1 • ;=1 /si 


M 


M 


M 


F l 2 =F 2l = Efa),fa)j> F,3=F 3 , =2(X,),.F s = F m =Efe)y. 

7=1 7=1 7=1 
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MM M 

G, =XM,S , , G 2 = £(*,), S, .“d G, = 2X 

7 = ■ 7=1 7=1 

The code developed (in MATLAB) for the present work has been given in Appendix B. In 
this, each signal is a set of 512 data points. At the end AjRj and A 2 R 2 give the separated 
modes of the signal. Figs. 4.5 and 4.6 show the signal modes, which are separated by 
least-squares method for the overlapped signals of Aluminium inserted Perspex sample 
and Epoxy resin inserted Perspex sample respectively. 

4.3 CROSS CORRELATION METHOD USING FFT 
In statistical terms, correlation is a measure of the linear relationship between two 
variables. The correlation between x and y can be written as, Zxy, where x and y are set 
of independent data points. It is also applicable to deterministic functions such as 
sinusoids. In the application of the correlation methods to engineering problems, the 
features generally expressed as waveforms rather than discrete time series. The 
correlation or covariance is expressed as a function of time delay r , as 

00 

C xy (r)= \x(t)y(t + r)dt (4.8) 

-00 

where, r is the time shift, lag or delay. 

If the waveform, x(t) correlate with a time-shift version of itself, is called 
Autocorrelation. Then the equation can be writing as 

00 

C x (t)= jx(t)x(t + r)dt (4.9) 

—00 

When a waveform correlate with another waveform, is called cross correlation, which 
describes the dependence between those waveforms. Then, the cross correlation function 
between waveforms, x(t) and y(t) can be expressed in terms of time shift as given in the 
equation 4.8. 

Correlation using FFT 

The correlation between the signals R(t) and S(t) can be written as 
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(4.10) 


C^(r)= \R(t).S{t + r).dt 
0 

where, R(t) is a reference signal, S(t) is the multi path signal and T is the time range. 
After Fast Fourier Transform (FFT), the equation 4.2 becomes 

S(v) = Yj A > ex P (-2-7 .v/) (4.11) 

/=1 

Then, the equation for p th mode can be written as 

F M = S(v)- - A, ,R(v). exp(- 2jjc.t,v) (4. 1 2) 

/=1 

where, s p (v) is the p th separated mode, S(v) and R(v) are the FFT of the S(t) and R(t). v is 
the frequency array, j indicates complex number. N is the number of modes. The 
separation of the various modes of the multi path signal has been done as per the flow 
chart as shown in fig. 4.7. At the end of the each iteration, the values of amplitude and 
time-shift are updating as per the equations 



and, 

z p = />os[max(F~' (s p (v)i?(v))) (4. 14) 

where, F 1 indicates Inverse Fast Fourier Transform (IFFT), max is for maximum value, 
pos for position for the maximum value index. These values of A p and x p will be using in 
the next iteration according to the flow chart [5] shown in fig. 4.7 for further updating. 

After finding the values of amplitudes and time-shifts for all the modes (N) of the 
signal, the individual modes can be reconstructed using amplitudes (A p ), time shifts (t P ) 
and reference signal (in frequency mode) as per the equation 

s, = R(v).A , . exp(-2 j.n xx t .v), i = 1,2,3, N (4. 1 5) 

where, N is the number of modes in the multi path signal. Appendix C shows the code 
written in MATLAB language for the present work. Figs 4.8 shows an example of the 
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reference signals, multi path (signal to be separated) signals, and their separated modes of 
the signals for sample of Aluminium inserted Perspex sample. In this, the overlapped 
signal is overlap of two signal modes with large time shift. Fig. 4.9 shows for overlapped 
signal with overlap of two signal modes with small time shift. Fig. 4.10 shows the 
separation of overlapped signal with very small time shift when the thick ness of the 
insertion is reduced. Figs 4.11 and 4.12 show the reference signals, multi path (signal to 
be separated) signals, and their separated modes of the signals for sample of Perspex with 
elliptical insertion of Perspex material. 

4.4 SELECTION OF DESIRED MODE 

After separation of individual modes of a multi-path signal, selection of desired mode 
from the separated modes of the signal plays major role in the Time of Flight (TOF) 
measurement. In the present work, a measure of energy criterion has been proposed for 
this selection. The measure of energy of a signal is defined as the sum of the squares of 
the signal values. Then the equation of the energy (EN) of the signal (s) having M data 
points is given as 

M 

i» = 2>, 2 (4.16) 

/ = 1 

In the present work, the measure of energy of the signals through base material 
and through insertion are differing a lot. Different attenuation characters of the material 
and scattering effects at the solid-solid interface may attribute this. So, for selection of 
desired mode from the separated modes, percentage of energy (normalized energy) of the 
signals have been taken into consideration rather than energy of the signals. Fig 4.16 and 
fig 4.17 show the energy of the signal versus probe position and percentage of energy 
versus probe position respectively. At position 15.5 mm the normalized (or percentage 
of) energy of the signal mode through insertion is more than the normalized (or 
percentage of) energy of the signal mode through base material. According to the energy 
criterion, the signal mode through insertion has been selected as a desired mode for the 
time of flight data. This energy-based criterion is equally applicable to the situations 
when the sound velocity in insert material is higher than the base material and vice versa. 
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4.5 TIME OF FLIGHT ESTIMATION 

Time of flight is the important parameter for reconstruction of the ultrasonic image with 
computer-aided tomography. The equation used for Time of flight data the equation as 
follows 

TOF = ptd + t (4.17) 

where, TOF is the Time of flight in second, ptd is the post trigger delay and t is the time 
delay (from an digitized signal). Fig 4.13 shows the post trigger delay and time delay 
measurements for TOF. In multi path signal, when the normalized energy of the 
separated signal mode received through insertion is higher than the signal through the 
base material, then the signal through insertion material has been taken into consideration 
for Time of flight estimation as shown in figs. 4.14 and 4.15 for the sample, sound 
velocity of insertion material is lower than base material and for the sample, sound 
velocity of insertion material is higher than base material respectively. 

4.6 CLOSURE 

In this Chapter about resolution of multi path signal using least squares method and cross 
correlation method using FFT are discussed. The energy criterion for selection of desired 
mode from the separated modes for TOF measurement is also presented. 
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Fig. 4.1(a) Sample- 1 of Perspex (B) with 20mmx20mm insertion of Epoxy Resin (A) with three 
positions of data collection 
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Fig. 4.1(b) Signals at three different positions (i.e. shown in fig. 4.1(a)) for Perspex sample with 
Epoxy resin insertion of 20mm x20mm 
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Fig. 4.1(c) Signals at three different positions (i.e. shown in fig. 4.1(a)) for Nylon sample 
with Epoxy resin insertion of 20mm x20mm 
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Fig. 4.1(d) Signals at three different positions (i.e. shown m fig. 4.1(a)) for Perspex 
sample with Aluminium insertion of 20mm x20mm 

























Fig. 4.2 (a) Sample-3 of Perspex (B) with 20mmxl0mm insertion of Epoxy Resin (A) 
with three positions of data collection 
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Fig. 4.2(b) Signals at three different positions (i.e. shown in fig. 4.2(a)) for Perspex 
sample with Epoxy resin insertion of 20mm xlOmm 
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Fig. 4.3(a) Sample-6 of Nylon (B) with 30mmx20mm elliptical insertion of Epoxy Resin 
(A) with four positions of data collection 
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Fig. 4.3(b) Signals at four different positions (i.e. shown in fig. 4.3(a)) for Nylon sample 
with Epoxy resin elliptical insertion 
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Fig. 4.3(c) Signals at four different positions (i.e. shown in fig. 4.3(a)) for Perspex sample 
with Epoxy resin elliptical insertion. 




Fig. 4.4(a) Sample-6 of Nylon (B) with 30mmx20mm elliptical insertion of Epoxy Resin 
1 1 (A) with four positions of data collection when rays passing through the 

(B) sample at 90° 
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Fig. 4.4(b) Signals at four different positions (i.e. shown in fig. 4.4(a)) for Nylon sample 
with Epoxy resin elliptical insertion. 
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Fig. 4.4(c) Signals at four different positions (i.e. shown in fig. 4.4(a)) for Perspex sample 
with Epoxy resin elliptical insertion. 
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Fig. 4.5 Separation of Overlapped Signal with two Reference signal using least square 
method for Perspex sample with Aluminium rectangular insertion (20mmx 
20mm) 



Fig. 4.6 Separation of Overlapped Signal with two Reference signal using least squares 

method for Perspex sample with Epoxy resin rectangular insertion (20mmx 
20mm) 
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Fig 4.7. 



Where, 


/ = // 


Separated signal mode, S p (v) = 5(v) — A t . exp(— 2 inx l v) 


/=i 

i*p 


Amplitude, A 


k(v)i 2 rfv 


, and 


I 

Time shift, r p = ^[max^FFT^ (v)i?(v 


Flow chart for the mode separation algorithm of a signal with linear 
superposition of temporarily shifted reference signal 
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Fig 4.8. Separation of Overlapped Signal with a Reference signal using cross correlation 
method for Perspex sample with Aluminium rectangular insertion (20mm x 
20mm) 



Fig 4.9. Separation of Overlapped Signal with a Reference signal using cross correlation 

method for Perspex sample with Epoxy resin rectangular insertion (20mm x 
20mm) 
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Fig 4.10. Separation of Overlapped Signal with a Reference signal using cross correlation 

method for Nylon sample with Epoxy resin rectangular insertion (20mmx 1 0mm) 



Fig 4. 1 1 . Separation of Overlapped Signal with a Reference signal using cross correlation 
method for Nylon sample with Epoxy resin elliptical insertion of 30mm major 
diameter and 20mm minor diameter, at insertion border 
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Fig 4. 12. Separation of Overlapped Signal with a Reference signal using cross correlation 
method for Nylon sample with Epoxy resin elliptical insertion of 30mm major 
diameter and 20mm minor diameter, at 2mm inside of insertion 



Fig 4. 13. Time-delay (t) and post trigger delay (ptd) measurements for the Time of Flight 
(TOF) data from the signal (Time v/s Amplitude) 
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Fig 4. 14. Time-delay (t) and post trigger delay (ptd) measurements for the Time of Flight 
(TOF) estimation from an overlapped signal, when wave velocity in insert 
material is slower than in base material. 
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Fig 4. 1 5 . Time-delay (t) and post trigger delay (ptd) measurements for the Time of Flight 

(TOF) estimation from an overlapped signal of sample, when wave velocity in 
insert material is faster than in base material. 
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Fig 4.16. Energy of signals through base-material and through insertion verses various 
positions from one of the edges to center of the insertion for Perspex sample with 
Epoxy resin rectangular insertion of 20mm x 20mm 



Fig 4. 17. Percentage of Energy (i.e. Normalized Energy) of signals through base- 
material and through insertion verses various positions from one of the edges 
to center of the insertion for Epoxy sample with Epoxy resin rectangular 
insertion of 20mm x 20mm 
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Chapter 5 

ULTRASONIC TOMOGRAPHIC RECONSTRUCTION 

5.1 INTRODUCTION 

Computerized Tomography (CT) is widely used in large number of scientific, medical 
and technical fields. In conventional projection radiography, the result is obtained just in 
the form of a shadowgarph requiring a good interpretation. CT, however, uses 
tomographic reconstruction to combine the information from multiple projections. Thus, 
the CT information is displayed quantitatively as a reconstructed slice plane of the object 
and it can provide accurate quantitative measures of material characteristics in small 
volume of the component. 

Although Computer Tomography was conceived originally for X-rays, in recent 
times, CT has been implemented successfully for other energy sources such as ultrasonic, 
laser, magnetic etc. Among these, Ultrasonic Tomography (UT) has special importance in 
application to Non Destructive Evaluation (NDE). Ultrasonic tomography can be defined 
as a technique to determine values of a spatially varying acoustic parameter of interest 
across a slice through an object. 

5.2 RECONSTRUCTION TECHNIQUES 

Reconstruction techniques are broadly classified into two groups: 

1 . Transformation methods or direct methods. 

2. Series expansion methods. 

Direct methods include direct algebraic methods and Fourier Transform methods and 
series expansion methods include Algebraic Reconstruction Technique (ART), 
Multiplicative Algebraic Reconstruction Technique (MART), Simultaneous Iterative 
Reconstruction Technique (SIRT), Entropy Optimization etc. ART and MART are 
further classified into different classes depending upon the way corrections are applied to 
the algorithms, are described below. 
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5.2.1 ART 

The iterative process starts from an initial approximation of slowness (reciprocal of the 
velocity) s° and the current iterate s k is corrected to a new iterate s k+l by taking into 
account only a single ray, say i‘ h , and changing only the image values of the pixels which 
intersect this ray. The discrepancy between the measurement of TOF, r, and the pseudo- 

N 

projection data Sj) obtained from the current image s is redistributed 

j = i 

among the pixels along the i ,h ray proportionally to their weights w y (i.e. intercept of the 

i ,h ray in the j th pixel) in the whole ray. In this way, the pixel values along the i' h ray are 
corrected to conform to the i‘ h measurement without changing the rest of the image. The 


whole process may be expressed as 


. . . T —< W' ,S k > 

s k+l =s k +- w‘ 

II W* II 2 

(5.1) 

where, 


<w i ,s k >=Y J w y sk j 

(5.2) 

N 

w 1 =Y J w IJ 

J = 1 

(5.3) 

and 


j|w'|| 2 =< w‘ ,w' > 

(5.4) 


Introducing the relaxation parameter, the modified step for updating the iterate can be 
expressed as 


s k+l = s k + A 


T~< w',5* > , 

— — — u; 

II w' II 2 


(5.5) 


where, A is the relaxation parameter, which is the real number, usually confined to the 
interval. 


£i <A<2 -e 2 , e x s 7 >0 


(5.6) 
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If the correction in negative, it may happen that the calculated field of a reconstructed 
image negative, in which case it is set to zero, as discussed by Gordon. 

ART may be further classified in different categories depending upon the correction 

applied. The four forms of ART are mentioned below. For all these algorithms, let ,v, for 

j=l,2...N be the initial approximation of field value or pixel value. Where N is the total 
no of pixels. 

Four classes of ART are 

1. MayingerART 

2. Gordon ART 

3. Gilbert ART 

4. Anderson ART 

5.2.2 MART 

When the correction is multiplicative, the ART is called Multiplicative AR T (MART). 
While in MART, at each iteration, a new estimate of field value is calculated by 
multiplication of previous value and a correction term, i.e. 

*r-Cx,; j = 1,2, N. (5.7) 

where C is correction term, j y is the slowness value of the /’ pixel and k is the iteration 

number. Depending up on the correction term, the MARI' may be further classified in 
four classes. Algorithm remains the same for all but correction term varies. Algorithm of 
MART may be stated as follows 
For each iteration k 

Step 1 : For each ray i , calculate the approximation TOP r, 

N 

T ' 

y-i 

Step 2: Modify the value of J* pixel, for all j* 1 ,2 ... N, as 

s* = Cxs* 1 


(5.K) 
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th 

Step 3: Repeat steps 1 and 2 for all rays. This completes the k iteration. 

Step 4: Repeat steps 1 to 3 until 


„i+l „k 

s — s 


x 100 < e 


(5.10) 


where e is the stopping criteria, generally take it as error. 


Four classes of MART are 

1) GBH MART 

2) GH MART 

3) Lent MART 

4) Lent2 MART 


5.2.2. 1 GBH MART 

Gordon, Bender, and Herman (GBH) first proposed a MART algorithm for CT 
applications. Like all MART algorithms, it is an iterative technique. In this method the 
correction term is given by 


C = 


RI 


V T :J 


if r, * 0, andw tJ * 0 
else 


(5.11) 


where, 

X is Relaxation parameter, 

* "I * ^ 

w tJ , is the Weight, intercept length in the j th pixel with i ray, 

{w lJ ) max , is the maximum of all weights on the i‘ h ray, 

r i ,r l , are measured and calculated projection data (TOF) respectively. 
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5.2.22 GH MART 

This MART algorithm was proposed by Gordon and Herman (GH), with a slight and 
important modification of the GBH algorithm. The correction term is given by, 

if r, * 0 ,andw v *0 (5.12) 

else 

5.2.23 Lent MART 

Lent has proposed another multiplicative updating scheme, in which the correction term 
is computed as, 

tfr,* 0, (5.13) 

else 

5.22A Lent2 MART 

A slight modification leads to Lent2 MART. The correction term is given by, 


if*,* 0, (5.14) 

else 


* rK„ 

where w tJ is given by — - 
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5.3 CONVERGENCE CRITERION 

It is necessary to determine when an iterative algorithm has converged to a solution, 
which is optimal. In all algorithms, iterations are carried out until the following relation is 
satisfied, 


k + 1 k 
S ~S 


xl00<e 


(5.15) 


where, e is the stopping criterion. In the present work, e has been taken as 0.001. Here 
/ +1 is the slowness (field value) at (k + 1)'* iteration and s k is the slowness at k th 
iteration. 

5.4 RESULTS AND DISCUSSION 

In the present work, Gordon and Herman (GH) MART algorithm has been used for 
reconstruction of CT images, which essentially computes the slowness (inverse of 
velocity) distribution in the domain from numerous projections (TOF) data from each 
different view. In this work, ultrasonic waves travel through the sample (from a series of 
source locations) to receivers (detectors) at the opposite positions on the other side of 
the object in the through transmission method. In the Ultrasonic Computer Tomography 
(UCT), two types of approaches are using. They are straight-ray approach and curved-ray 
approach. In straight-ray approach the ray paths are assumed as straight lines. But, in the 
curved-ray approach ray-paths are assumed, with little bending. In the present work, 
simple straight-ray approach has been used. 

The present study has been conducted on domains of various samples containing 
different inclusions as given in the following Table. 5.1. 


gwfrrc #sn?r 

a.. 



»«+***#*• IM»« 
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Table 5.1 


Sample 

Number 

Material and Dimensions (All in mm) 

Longitudinal wave 
Velocity (m/s) 

Projection 

datatype 

Base 

Insert 

Body 

Insert 

1 . 

Perspex 

(50x50 Square) 

Epoxy Resin 

(30x20 Ellipse) 

2730 

2650 

4-view 

2. 

Perspex 

(50x50 Square) 

Epoxy Resin 

(20x20 Square) 

2730 

2650 

4- view 

3. 

Nylon 

(50x50 Square) 

Epoxy Resin 

(30x20 Ellipse) 

2300 

2650 

2-view 


For reconstruction purpose the domain of each sample is divided into 50x50 
grids (pixels, i.e. Cartesian grid of square picture elements) as shown in fig. 5.1(b). The 
pixels were numbered from 1 to 2500 (i.e. 50 x 50). While ultrasonic signal is passing 
through the material, the ray crosses various pixels. When the i‘ h ray is passing through 
the j th pixel, the length of the intercept of the i ,h ray in the j th pixel is denoted by w tj , 

called weight, as shown in fig. 5.1(a). In each view (of angle), the TOF data has been 
taken with 1mm of interval, called ray spacing. In 4- view projection, the analysis was 
repeated for four views (directions), namely 0°, 90°, 45° and 135°, with respect to positive 
x-axis as shown in figs. 5.1(c) through 5.1(f). But, in case of 2-view reconstruction, the 
data has been considered for 0° and 90° views only. In the present work, 50 rays have 
been considered for both 0° and 90°, and in each view, 71 rays have been considered for 
both 45° and 135°. In the present work, for getting the data for 45° and 135°, some special 
setup was used as shown in fig. 2.2(h). The ray locations were so chosen that the ray 
spacing (perpendicular distance between rays) is uniform for all views (i.e. one pixel 
dimension, 1mm). TOF data has been computed for all rays, as discussed in the previous 
chapter. 
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The reconstruction was conducted from projection (TOF) data. Projections have 
been taken either from four views (directions) namely 0°, 90°, 45° and 135° (in case of 4- 
view projection) or from two views (directions) namely 0° and 90° (in case of 2- view 
projection). The details of each view are discussed below. 

1. Ray-paths along 0° direction (x-axis) 

In this view, in all 50 data have been collected from 50 ray-paths, each ray originating 
from a pixel in the first column (fig. 5.1(c)) of the grid and being excited parallel to x- 
axis. In this case first ray-path has been taken at 0.5mm (approximately) distance from x- 
axis and parallel to x-axis. The other rays have been passed parallel to the x-axis with 
uniform ray spacing, being equal to unit pixel length (i.e. 1mm). The corresponding 
receiver to any receiver to any source, in this case, was situated exactly opposite it in the 
x-direction, as shown in fig. 5.1(c). 

2. Ray-paths along 90° direction (x-axis) 

As explained earlier, in this view also 50ray-paths have been taken into consideration. In 
this case, each ray originating from every column of the grid, towards y-direction as 
shown in fig. 5.1(d). In this case also the ray spacing (i.e. 1mm) is constant and is same 
as that in case of 0° view. 

3. Ray-paths along 45° direction 

In this case, 71 ray-paths have been considered. In this case also the ray spacing is 
uniform and is equal to 1mm. But, for getting TOF data, some special setup has been 
used as shown in fig. 2.2(h), by using four right-angle triangular Perspex (same as base 
material of the specimen) pieces. The time of flight was taken for this extra material, has 
been subtracted from the experimental TOF data for the getting input projection (TOF) 
data. 

4. Ray-paths along 135° direction 

As before, 71 ray-paths have been considered in this view as shown in fig. 5.1(f), with 
ray spacing of 1mm (same as above cases). 
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In all these cases, some water gaps (approximately 1mm each) were maintained 
between source probe and specimen plate surface, and between specimen plate surface 
and detector probe, for free movement of the probes. The time of flight for the signal 
through these gaps were taken into recount in analysis the experimental TOF data. 

5.5 IMAGE RECONSTRUCTION 

In the present work of Ultrasonic Computer Tomography (UCT), a total of 242 projection 
data are available for 4- view reconstruction, where as 100 data are available for 2-view 
reconstruction. The reconstruction algorithm here determines the slowness value of each 
pixel. 

The relation between the TOF, slowness and weight (intercept length) matrix is 
given by. 

[(F]{S}={r} (5.16) 

where, W is the weight matrix, of 242 ray x 2500 pixels in case of 4-view reconstruction 

and 100 ray x 2500 pixels in case of 2-view reconstruction, S is slowness vector of 2500 
pixels, and t is the TOF vector of 242 ray-paths in case of 4-view and of 100 ray-paths in 
case of 2-view reconstruction. 

The reconstruction procedure has been implemented in four steps. 

1 . In this first step, the algorithm determines the co-ordinates of the all pixels of the 
domain. 

2. In this step the algorithm evaluates the weight matrix of 242 x 2500 (or) 100 x 2500 
dimension, with help of first step results. 

3. In this step, using weight matrix and projection data (Time of flight data), it finds the 
slowness vector using GH MART algorithm. In this program, the first iteration starts with 
initial guess pixel (slowness) values, other than zeros as, it is multiplicative algorithm. In 
this reconstruction algorithm, the pixel values are updates to bring in to better agreement 
with the measured projection. This updating of pixel values is continued until a stopping 
criterion is met, as shown in equation.5.10. A possible stopping criteria is the absolute 
incremental and the change is less than some prescribed value, the iteration process 
stopped. 
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4. Finally, the slowness values all pixels were changed to the values between 0 and 255 
in PGM format to plot the reconstructed image. 

In the present work, Normalized RMS Error (L2 error) has been used to find the 
deviation of reconstructed field from the exact field. The Normalized Root Mean Square 
Error (EL 2 ) is defined as 



V — V 
max mm 


(5.17) 


where, s l a , s l r are the actual and reconstructed slowness values of the i lh pixel 
respectively. N is total number of pixels. s miX , s min are maximum and minimum of all 
pixel values respectively. The actual slowness of the Perspex, Nylon and Epoxy resin 
have been taken as 366.3 jus/m 2 , 434.8 juslm 2 and 381.7 jjs/m 2 respectively. 

For all the samples, the convergence criterion has been taken as 0.001. The 
percentage of sound velocity difference between the base material and the insert material 
is less than 16% in all the samples, so straight-ray assumption is valid. In the present 
work, the effects of the relaxation parameter on the reconstruction accuracy and 
convergence rate have also been presented. 


(i). Sample 1: 

The sample is made of Perspex with rectangular insertion of Epoxy resin and has been 
tested with different relaxation parameters. Table 5.2 shows the Normalized RMS (EL 2 ) 
error and number of iterations till the convergence criterion is met at various relaxation 
parameters (X). From these results, it has been found that, at relaxation parameter X=0.35, 
it gives better accuracy of the reconstructed image (i.e. EL 2 error is minimum, 18.753%), 
with least number of iterations (i.e. 34). Datta [23] reported that 47.1% RMS error was 
the least at A,=1.0 for rectangular insertion sample with 4- view projection with FEM 
simulated data. Figs. 5.2(a) to 5.2(d) show the reconstructed image at four different 
relaxation parameters i.e. A.O.02, 0.35, 0.5 and 1.8. The insert is also shown by a white 
line. 


55 



ii. Sample 2: 

The sample is of Perspex as a base material with elliptical insertion of Epoxy resin. In 
this sample, the velocity of sound in insert material is slower than the base material. The 
accuracy of the reconstructed images in terms of Normalized RMS error (EL 2 ) and 
convergence rate has been tested at all different relaxation parameters (A.). Table. 5.3 
shows the EL 2 error and the number of iterations at various relaxation parameters. It can 
be seen that at A=0.316, the reconstructed image has minimum Normalized RMS error 
(i.e. 24.6722%) and least number of iterations (36). Figs. 5.3(a) to 5.3(d) show the 
reconstructed images at four different relaxation parameters, namely A.=0.05, 0.3, 0.316 
and 1.5. It shows that the image for A,=0.316, has better appearance compared to the 
reconstructed images at other relaxation parameters. It has also been found that, the 
Normalized RMS error depends on the shape and the size of the defect, as the EL 2 error 
of the elliptical inserted sample is more than the rectangular inserted sample, even though 
both samples were tested with 4- view projection data. 

iii. Sample 3: 

In this case, the Nylon sample with elliptical Epoxy resin insertion has been tested with 
2-view projection data, i.e. 0° and 90°. The accuracies of the reconstructed images and 
convergence rates have been tested at various relaxation parameters. Table 5.4 shows the 
Normalized RMS (EL 2 ) error and iterations at different relaxation parameters. In this test, 
it has been found that, the normalized RMS error (i.e. 32.4328) is minimum at A,=T.O with 
minimum number of iterations (i.e. 2). Figs. 5.4(a) to 5.4(d) show the reconstructed 
images at four relaxation parameters, i.e. A,=0.01, 0.1, 1.0 and 1.8. 

From above tests for the three samples using the present method of TOF estimation the 
reconstruction seem to be better. It can also seen that, reconstruction from 4-view 
projection data gives the better results than 2- view reconstruction. And it has also been 
found that, at optimum relaxation parameter, the reconstruction gives more accurate 
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Fig. 5.1(c) Ray paths along 0° direction (x-direction) 

Y 



SOURCES 

Fig. 5.1(d) Ray paths along 90° direction (y-direction) 




(a) X = 0.02 



(b) X = 0.35 




(c) X = 0.5 


Figs. 5.2(a-d) 



(d) X = 1-8 


from 4- view projection data for 
relaxation parameters (A.) for Epoxy 
(faster material) sample. 


Tomographic reconstructed images 
Rectangular Insertion with different 
(slower material) inserted in the Perspex 
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(c) A. = 0.316 (d) A,= 1.5 


Figs. 5.3 (a-d). Tomographic reconstructed images from 4-view projection data for 
Elliptical Insertion with different relaxation parameters (1) for Epoxy 
(slower material) inserted in the Perspex (faster material) sample. 


(a) X = 0.05 


(b) X = 0.3 
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(c) A, = 1 .0 (d) X = 1.8 

Figs. 5.4 (a-d). Tomographic reconstructed images from 2-view projection data for 
Elliptical Insertion with different relaxation parameters (X.) for Epoxy 
(faster material) inserted in the Perspex (slower material) sample 
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s. 

No. 

Relaxation 
Parameter ( [X ) 

Number of 
Iterations 

Normalized 
RMS Error 
(EL2) 

1 

0.02 

257 

18.853 

2 

0.05 

136 

18.853 

3 

0.08 

96 

18.853 

4 

0.10 

82 

18.847 

5 

0.20 

50 

18.818 

6 

0.30 

38 

18.774 

7 

0.34 

35 

18.759 

8 

0.345 

35 

18.755 

9 

0.348 

35 

18.754 

10 

0.349 

35 

18.754 

11 

0.35 

34 

18.753 

12 

0.36 

34 

18.755 

13 

0.4 

252 

19.281 

14 

0.5 

422 

20.023 

15 

1.0 

498 

24.393 

16 

1.5 

318 

25.101 

17 

1.8 

197 

24.815 


Table 5.2 Normalized RMS Error of the reconstructed images of die Rectangular 
Epoxy inserted Perspex sample from 4-view projection data with different 
relaxation parameters (X) 
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S.No. 

Relaxation 
Parameter (A,) 

Number of 
Iterations 

Normalized 

RMS Error (EL2) 

1 

0.05 

125 

24.87 

2 

0.08 

88 

24.85 

3 

0.1 

• 74 

24.8377 

4 

0.2 

47 

24.7746 

5 

0.3 

37 

24.6884 

6 

0.31 

36 

24.6792 

7 

0.311 

36 

24.6765 

8 

0.312 

36 

24.6757 

9 

0.313 

36 

24.6745 

10 

0.314 

36 

24.6741 

11 

0.315 

36 

24.6729 

12 

0.3155 

36 

24.6724 

13 

0.316 

36 

24.6722 

14 

0.317 

36 

24.6734 

15 

0.32 

157 

24.8908 

16 

0.35 

283 

25.1578 

17 

0.4 

414 

25.5147 

18 

0.5 

534 

26.05 

19 

1.0 

541 

27.41 

20 

1.5 

355 

28.248 


Table 5.3 Normalized RMS Error of the reconstructed images of the Epoxy elliptical 
inserted Perspex sample from 4-view projection data with different 
relaxation parameters (IV) 
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S. No. 

Relaxation 
Parameter (X) 

Number of 
Iterations 

Normalized 

RMS Error (EL2) 

1 


438 

32.407 

2 

0.02 

254 

32.419 

3 

0.05 

119 

32.423 

4 

0.08 

79 

32.423 

5 

0.1 

65 

32.423 

6 

0.5 

14 

32.433 

7 

0.8 

7 

32.4333 

8 

0.9 

5 

32.4333 

9 

0.95 

5 

32.4328 

10 

1.0 

2 

32.4328 

11 

1.05 

5 

32.4328 

12 

1.1 

5 

32.4328 

13 

1.2 

7 

32.4333 

14 

1.5 

15 

32.4333 

15 

1.8 

46 

32.42226 


Table 5.4 Normalized RMS Error of the reconstructed images of the Epoxy elliptical 
inserted Nylon sample from 2- view projection data with different relaxation 
parameters (X) 
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(a) A, = 0.1 


(b) X = 0.5 



Fig. 5.5 (a-b) Tomographic reconstructed images from 4-view projection data for. 

Rectangular Insertion with different relaxation parameters (A.) for Epoxy 
(slower material) inserted in the Perspex (faster material) sample with 
out signal mode separation 


S. No. 

Relaxation 
Parameter (A.) 

Number of 
Iterations 

Normalized 

RMS Error (EL2) 

1 

0.1 

92 

21.307 

2 

0.5 

538 

36.842 

3 

1.0 

967 

38.644 

4 

1.5 

1023 

38.851 


Table 5.5 Normalized RMS Error of the reconstructed images of the Rectangular 
Epoxy inserted Perspex sample from 4-view projection data with different 
relaxation parameters (A.) with out signal mode separation 
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(a) A, = 0.5 (b) A, = 1.5 

Fig. 5.6 (a-b) Tomographic reconstructed images from 4-view projection data for 
Elliptical Insertion with different relaxation parameters (A.) for Epoxy 
(slower material) inserted in the Perspex (faster material) sample with 
out signal mode separation. 


S. No. 

Relaxation 
Parameter (A) 

Number of 
Iterations 

Normalized 

RMS Error (EL2) 

1 

0.1 

76 

25.731 

2 

0.5 

617 

26.529 

3 

1.0 

608 

29.492 

4 

1.5 

838 

32.711 


Table 5.6 Normalized RMS Error of the reconstructed images of the Epoxy elliptical 
inserted Perspex sample from 4-view projection data with different 
relaxation parameters (A) with signal mode separation. 
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(a) A, = 0.1 


(b) A, = 1 .0 


Fig. 5.7 (a-b) Tomographic reconstructed images from 2-view projection data for 
Elliptical Insertion with different relaxation parameters (A.) for Epoxy 
(faster material) inserted in the Perspex (slower material) sample with 
signal mode separation. 


S. No. 

Relaxation 
Parameter (A.) 

Number of 
Iterations 

Normalized 

RMS Error (EL2) 

1 

0.1 

65 

33.183 

2 

0.5 

14 

33.151 

3 

1.0 

2 

33.121 

4 

1.5 

15 

33.161 


Table 5.7 Normalized RMS Error of the reconstructed images of the Epoxy elliptical 

inserted Nylon sample from 2-view projection data with different relaxation 
parameters (A,) with out signal mode separation. 
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Chapter 6 

CONCLUSIONS AND SCOPE FOR FUTURE WORK 

6.1 CONCLUSIONS 

In the present work, ultrasonic tomographic reconstruction is used for reconstruction of 

domains of the samples with artificially implanted insertions. Based on the results 

obtained, the following conclusions are drawn: 

• The wavelet transform is an excellent tool for noise suppression. A feature of this is 
that noise can be suppressed by utilization of different levels; where each level is 
depended upon the type of the signal. 

• Correlation technique can be effectively employed for mode separation of the 
overlapped signal with any number of modes. This analysis results in accurate 
estimation of time of flight. 

• A measure of energy-based criterion is proposed for selection of desired mode from 
the separated modes of an overlapped signal. This criterion is equally applicable for 
both fast and slow insertion medium. 

6.2 SCOPE FOR THE FUTURE WORK 

• These noise suppression and signal separation techniques can be extended to defect 
characterization in composite materials. 

• More studies have to be made on wavelet analysis to apply it for mode separation in 
the overlapped signals. 

• By the reconstruction of the images for 2-D slices of the defective sample, the 
reconstruction of 3-D defects can be carried out by integrating the 2-D slice images. 
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Avnendix A 


PROGRAM FOR NOISE SUPPRESSION 

Wavelet toolbox is a collection of functions on the MATLAB technical computing 
environment. It has been used for noise suppression of a given signal using ‘db5’ 
wavelet. The program used for the present work is given below. 

% MATLAB PROGRAM FOR WAVELET DECOMPOSITION FOR NOISE 
% SUPPRESION OF A SIGNAL WHICH IS USED IN THE PRESENT WORK 

wname- db5'; % Name of the wavelet which is using 

NUM=4; % Number of the signal for decompose 

% Reading the data from a file 
format short e; 

fid=fopen('FILENAME.m7r'); 
n=fscanf(fid,'%d',[l .1]); 

AA=zeros(n,5 1 2);B=zeros( 1,512); 
for i=l:n 

[B,count]=fscanf(fid,'%d',[l,512]); 

AA([i],:)=AA([i],:)+B; 

end 

status=fclose('aH') ; 

AA; 

A=AA([NUM],:); 

plot(A); 

pause 

% Wavelet decomposition using ‘db3’ wavelet 3 -level decomposition 

[C,L]=wavedec(A,3, wname); % Multi(3)-level 1-D wavelet decomposition 

cA3=appcoef(C,L,wname,3); % cA3 -Extract 1-D Approximation coefficients 

cD3=detcoef(C,L,3); % cD3-Extracts the Detail coefficients at level 3 

cD2=detcoef(C,L,2); % cD2-Extracts the Detail coefficients at level 2 

cD 1 =detcoef(C,L, 1 ); % cD 1 -Extracts the Detail coefficients at level 1 

% Reconstruct from Approximation coefficients-cA3 
A3=wrcoef('a',C,L,wname,3); 

% Reconstruct from Detail coefficients-cD3 
D3=wrcoef('d',C,L,wname,3); 

% Reconstruct from Detail coefficients-cD2 
D2=wrcoef('d',C,L,wname,2); 
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% Reconstruct from Detail coefficients-cDl 
D 1 =wrcoef('d',C,L 5 wname, 1 ) ; 

% Plotting the signal Approximation and details 
subplot(5,l,l);plot(A);gtext('Original Signal S') 
subplot(5, 1 ,2);plot(A3);gtext(' Approximation A3'); 
subplot(5, 1 ,3);plot(D3);gtext('Detail D3 ') 
subplot(5,l ,4);plot(D2);gtext('Detail D2') 
subplot(5, 1 ,5);plot(D 1 );gtext('Detail D 1 ') 

% In this program ‘A3’ is the Signal after suppress the noise 
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Ajo^endbc B 


PROGRAM FOR SIGNAL SEPARATION USING LEAST-SQUARES 
METHOD 

In this program Least squares method used for separating overlapped signal in to two 
modes, when two-reference signal are exist. MATLAB environment has been used in the 
present work. The program is given below. 

<y Q ** He******************************************************** ********** 

% PROGRAM FOR SEPARATION OF THE OVERLAPPED SIGNAL USING 
% LEAST-SQUARES METHOD 

wname='db5'; % Name of the wavelet (i.e. using) 
fwavno=l; % First Reference signal index 
owavno=2; % Overlapped signal index (from input file) 
lwavno=3; % Second reference signal 

% Reading the data from a file 
format short e; 
fid=fopen('inputfile.m , ,'r'); 
n=fscanf(fid,'%d , ,[ 1 . 1 ]); 

DATA=zeros(n,5 1 2);B=zeros(l ,512); 


for i=l:n 

[B,count]==fscanf(fid,'%d',[l ,5 1 2]); 
DATA([i],:)=DATA([i],:)+B; 
end 

status=fclose('aH'); 

DATA; 

% Shifting the axis of the signal to the x-axis 
A=zeros(n,512); 
for i=l:n 

A([i] ,:)=DATA([i] , :)-round(mean(DATA([i] ,:))); 
end 


% Least squares method starts 

FT=A([owavno],:);Fl=A([fwavno],:);F2=A([lwavno],:); 

F 1 1 =sum(F 1 . * F 1 ) ;F i 2=sum(F 1 . * F2) ;F 1 3 =sum(F 1 ) ; 
F22=sum(F2.*F2);F23=sum(F2);F33=512; 

G 1 =sum(FT. * F 1 ); G2=sum(FT. * F2) ;G3 =sum(FT) ; 


73 



F=[F11 F12 F13;F12 F22 F23;F13 F23 F33] 


G=[G1;G2;G3] 

AA=F\G 

a=AA([l]);b=AA([2]);c=AA([3]); 
wav 1 =a* F 1 ; wav2=b*F2; 
twav=wav 1 +wav2+c ; 

% Plotting the signals 

subplot(5 ,1,1) ;plot(FT) ;gtext('Overlapped Signal') 
subplot(5,l,2);plot(Fl);gtext('First Reference Signal') 
subplot(5,l,3);plot(F2);gtext('Second Reference Signal') 
subplot(5,l,4);plot(wavl);gtext('First separated Signal') 
subplot(5,l,5);plot(wav2);gtext('Second separated Signal') 

% Finding the error of summation (of the separated) signal 
twav=wav 1 +wav2; 
e=0;T=0; 
for i=l:512 

al=FT([i]);bl=twav([i]); 
e=e+abs((al -b 1 )/b 1 ) ; 
end 


pause 

x=l :5 12;plot(x,twav,x,FT); 



Ajgyendix C 


PROGRAM FOR SIGNAL SEPARATION USING CROSS 
CORRELATION METHOD 

In this program Cross-correlation method used for separating overlapped signal in to two 
modes, when only one reference signal is exists. MATLAB environment has been used in 
the present work. The program is given below. 

(^ ********************************************************************* 
% PROGRAM FOR SEPARATION OF THE OVERLAPPED SIGNAL USING 
% CROSS CORRELATION METHOD 

0 /^ ********************************************************************* 

wname='db5'; % Name of the wavelet (i.e. using) 

fwavno=l ; % Reference signal index (in the input file) 

owavno=2; % Overlapped signal index 

% Reading the data from a file 
format short e; 
fid=fopen('inputfile.mVr'); 
n=fscanf(fid,'%d',[l . 1 ]); 

DATA=zeros(n,5 12);B=zeros(l ,5 12); 
for v=l:n 

[B , count] =fscanf(fid,'%d' , [ 1 ,512]); 

DATA([v],:)=DATA([v],:)+B; 

end 

status=fclose('aH') ; 

% Suppression noise of the signal 

DATAl=zeros(2,512); A3=zeros(2,512); 

DATAl=DATA([fwavno],:); DATA1 =DATA([owavno],:); 


for i=l:2 

[C,L]^avedec(DATAl([i],:),3,wname); 

c A3 =appcoef(C,L, wname, 3 ) ; 
cD3=detcoef(C,L,3); 
cD2=detcoef(C,L,2) ; 
cD 1 =detcoef(C,L, 1 ); 

A3 ( [i] , :)=wrcoef('a',C,L, wname,3 ); 
D3=wrcoef('d',C,L,wname,3); 
D2=wrcoef('d',C,L,wname,2) ; 

D 1 =wrcoef('d',C,L, wname, 1 ); 
end 

% Cross correlation program Starts from here 
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RT=A3([1],:);ST=A3([2],:); ^*^$5 

S=fft(ST,512); % FFT of the Overlapped signal 

R=fft(RT,512); % FFT of the Reference signa 

n=length(R); % size of the array 

RNU=R; SNU=S ; 
nyquest=l/2; 

NU=(0:n-l)/(n/2)*nyquest; % Frequency scale 


N=2; % N -> Number of Modes 

A=zeros(N, 1 );T AU=zeros(N, 1 ) ; 

for k=l :N % k -> For k modes, k = 1 to N 

for p=k:-l : 1 % p -> For on Order p mode, p = k to 1 

TEMP=0+j*0; 
for g=l:N 
if g~=p 

TEMP=TEMP+A(g)*exp(-2*j*pi*TAU(g)*NU); 

end 

end 

SP=SNU-RNU.*TEMP; 

ASP=abs(SP); 

RNU_R=RNU(5 1 2 1 : 1 ); 

[MAX,pos] = max(abs(ifft(SP * RNU_R))); 

TAU(p)=pos; 

ARNU=abs(RNU); 

A(p)=sum(ASP. A 2)/sum(ARNU. A 2); A(p)=sqrt(A(p)) 
end 

end 

A,TAU 


% Finding separated signal(in frequency mode) from amplitudes & time shifts 
SPl=RNU.*(A(l)*exp(-2*j*pi*TAU(l)*NU)); 
SP2=RNU.*(A(2)*exp(-2*j*pi*TAU(2)*NU)); 

subplot(4,l,l);plot((ST),'k');grid on;gtext('Overlapped Signal S'); 
subplot(4,l,2);plot((RT),'k');grid on;gtext('Reference Signal R'); 

% Changing from frequency-domain to time-domain 
SPTl=ifft(SPl);SPT2=ifft(SP2); 

subplot(4,l,3);plot(real(SPTl),'k');grid on;gtext('First Separated Signal S_l'); 
subplot(4, 1 ,4);plot(real(SPT2),'k');grid on;gtext('Second Separated Signal S_2'); 

% Addition of of two modes of the signal 
SPT12=SPT1+SPT2; . 

Pause; subplot( 1,1,1 );plot(real(SPT 1 2)) 
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